The goal for the first project is to compute the magnetostatic problem for a "real" world problem.
To this end, we will work with Alessio & Peter's motor geometry of a permanent magnet machine.
For now, let $\Omega$ denote the computational domain.
Then Maxwell's equations for magnetostatics are described by:
\begin{alignat}{2} \operatorname{curl} h &= j, \qquad && \text{in } \Omega, \label{eq:1}\\ \operatorname{div} b &= 0, \qquad && \text{in } \Omega, \label{eq:2}\\ n \cdot b &= 0, \qquad && \text{on } \partial\Omega, \label{eq:3} \end{alignat}On simply connected domains, the second equation implies the existence of a vector potential $a$ such that
\begin{align} \operatorname{curl} a &= b\label{eq:4} \end{align}So far, these equations are independent of the materials.
To incorporate material properties, we need to introduce constitutive equations, which will be different depending on the material type, such as air, iron, or magnet.
For now, assume that the magnetic field $h$ and magnetic flux $b$ are linked by
\begin{align} h = f'(b) \label{eq:nl} \end{align}where $f$ denotes a smooth and strongly convex "energy density" $f: \mathbb{R}^d \to \mathbb{R}$.
Using this relation, we can rewrite the problem in a compact manner in the following way:
\begin{align*} \operatorname{curl}(f'(\operatorname{curl} a)) &= j + \operatorname{curl} m \\ n\cdot \operatorname{curl} a &=0 \end{align*}In 2D, we can make the following simplifications:
In the typical setting with cylindrical symmetry, which is of relevance for electric motors, one has $\Omega = \Omega_{2D} \times (0,L)$ and assumes that $j=(0,0,j_3)$, $m=(m_1,m_2,0)$, $a=(0,0,a_3)$, $b=(b_1,b_2,0)$, and $h=(h_1,h_2,0)$, with all components independent of $z$.
In this way, one arrives at the scalar problem
\begin{alignat}{2} \operatorname{Curl} f'(\operatorname{curl}_{2D} a_3) &= j_3 + \operatorname{Curl}_{2D} m, \qquad && \text{in } \Omega_{2D}, \\ n \cdot \operatorname{curl}_{2D} a_3 &= 0, \qquad && \text{on } \partial\Omega_{2D}. \end{alignat}Here $\operatorname{Curl}_{2D}=(-\partial_y, \partial_x)$ and $\operatorname{curl}_{2D}=\binom{\partial_y}{-\partial_x}$ are the vector-to-scalar and scalar-to-vector curl, respectively.
This problem can further be transformed into
\begin{alignat}{2} -\operatorname{div} (g'(\nabla u)) &= j_3 + \operatorname{div} m^\perp :=\widehat{j_3}, \qquad &&\text{in } \Omega_{2D} \label{eq:sys1} \\ u &= 0, \qquad && \text{on } \partial\Omega_{2D} \label{eq:sys2} \end{alignat}where now $u=a_3$ is the $z$-component of the vector potential, $\binom{-b_2}{b_1}=\nabla u$ is the rotated magnetic flux, $\binom{-h_2}{h_1}=g'(\nabla u)$ is the rotated magnetic field, and $m^\perp = \binom{-m_2}{m_1}$ is the rotated magnetization.
The boundary condition for the transformed problem reads $\partial_t u := t\cdot\nabla u=0$, meaning that the tangential component of $u$ vanishes on the whole boundary. This, in turn, means that $u$ has to be constant on the boundary. For this reason, we may set $u=0$ on $\partial\Omega$. Careful! If the domain is not simply connected, we are not allowed to prescribe $u$ everywhere on the boundary!
Loading imports ...
from IPython.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
import sys; sys.path.insert(0,'../../') # adds parent directory
import time
import numpy as np
import pde
import scipy.sparse as sps
from sksparse.cholmod import cholesky as chol
import plotly.io as pio
pio.renderers.default = 'iframe'
import pyamg
import plotly
plotly.offline.init_notebook_mode()
ORDER = 1
Loading the motor, generating inital mesh
motor_npz = np.load('meshes/motor.npz', allow_pickle = True)
p = motor_npz['p'].T
e = motor_npz['e'].T
t = motor_npz['t'].T
q = np.empty(0)
regions_2d = motor_npz['regions_2d']
regions_1d = motor_npz['regions_1d']
m = motor_npz['m']
j3 = motor_npz['j3']
MESH = pde.mesh(p,e,t,q)
# MESH.refinemesh()
# MESH.refinemesh()
Prepares lists of indices.
def getIndices(liste, name, exact = 0, return_index = False):
if exact == 0:
ind = np.flatnonzero(np.core.defchararray.find(list(liste),name)!=-1)
else:
ind = [i for i, x in enumerate(list(liste)) if x == name][0]
elem = np.where(np.isin(MESH.t[:,3],ind))[0]
mask = np.zeros(MESH.nt); mask[elem] = 1
if return_index:
return ind, mask
else:
return mask
mask_air_all = getIndices(regions_2d, 'air')
mask_stator_rotor_and_shaft = getIndices(regions_2d, 'iron')
mask_magnet = getIndices(regions_2d, 'magnet')
mask_coil = getIndices(regions_2d, 'coil')
mask_shaft = getIndices(regions_2d, 'shaft')
mask_linear = mask_air_all + mask_magnet + mask_shaft + mask_coil
mask_nonlinear = mask_stator_rotor_and_shaft - mask_shaft
ind_stator_outer = np.flatnonzero(np.core.defchararray.find(list(regions_1d),'stator_outer')!=-1)
edges_stator_outer = np.where(np.isin(e[:,2],ind_stator_outer))[0]
ind_trig_coils = {}
for i in range(48):
ind_trig_coils[i] = getIndices(regions_2d, 'coil' + str(i+1), exact = 1, return_index = True)[0]
ind_trig_magnets = {}
for i in range(16):
ind_trig_magnets[i] = getIndices(regions_2d, 'magnet' + str(i+1), exact = 1, return_index = True)[0]
Setting up the function $f$ and its derivatives
k1 = 49.4; k2 = 1.46; k3 = 520.6; nu0 = 10**7/(4*np.pi)
# k1 = 3.8; k2 = 2.17; k3 = 396.2
f_iron = lambda x,y : k1/(2*k2)*(np.exp(k2*(x**2+y**2))-1) + 1/2*k3*(x**2+y**2) # magnetic energy density in iron
nu = lambda x,y : k1*np.exp(k2*(x**2+y**2))+k3
nux = lambda x,y : 2*x*k1*k2*np.exp(k2*(x**2+y**2))
nuy = lambda x,y : 2*y*k1*k2*np.exp(k2*(x**2+y**2))
fx_iron = lambda x,y : nu(x,y)*x
fy_iron = lambda x,y : nu(x,y)*y
fxx_iron = lambda x,y : nu(x,y) + x*nux(x,y)
fxy_iron = lambda x,y : x*nuy(x,y)
fyx_iron = lambda x,y : y*nux(x,y)
fyy_iron = lambda x,y : nu(x,y) + y*nuy(x,y)
f_linear = lambda x,y : 1/2*nu0*(x**2+y**2)
fx_linear = lambda x,y : nu0*x
fy_linear = lambda x,y : nu0*y
fxx_linear = lambda x,y : nu0 + 0*x
fxy_linear = lambda x,y : x*0
fyx_linear = lambda x,y : y*0
fyy_linear = lambda x,y : nu0 + 0*y
f = lambda ux,uy : f_linear(ux,uy)*new_mask_linear + f_iron(ux,uy)*new_mask_nonlinear
fx = lambda ux,uy : fx_linear(ux,uy)*new_mask_linear + fx_iron(ux,uy)*new_mask_nonlinear
fy = lambda ux,uy : fy_linear(ux,uy)*new_mask_linear + fy_iron(ux,uy)*new_mask_nonlinear
fxx = lambda ux,uy : fxx_linear(ux,uy)*new_mask_linear + fxx_iron(ux,uy)*new_mask_nonlinear
fxy = lambda ux,uy : fxy_linear(ux,uy)*new_mask_linear + fxy_iron(ux,uy)*new_mask_nonlinear
fyx = lambda ux,uy : fyx_linear(ux,uy)*new_mask_linear + fyx_iron(ux,uy)*new_mask_nonlinear
fyy = lambda ux,uy : fyy_linear(ux,uy)*new_mask_linear + fyy_iron(ux,uy)*new_mask_nonlinear
Setting first or second order variables, expanding mask for second order, as $f$ has to be evaluated for $P_1$, where the gradiant of our solution lies
if ORDER == 1:
poly = 'P1'
dxpoly = 'P0'
order_phiphi = 2
order_dphidphi = 0
new_mask_linear = mask_linear
new_mask_nonlinear = mask_nonlinear
if ORDER == 2:
poly = 'P2'
dxpoly = 'P1'
order_phiphi = 4
order_dphidphi = 2
new_mask_linear = np.tile(mask_linear,(3,1)).T.flatten()
new_mask_nonlinear = np.tile(mask_nonlinear,(3,1)).T.flatten()
Assembling the matrix and right-hand sides
tm = time.monotonic()
phi_H1 = pde.h1.assemble(MESH, space = poly, matrix = 'M', order = order_phiphi)
dphix_H1, dphiy_H1 = pde.h1.assemble(MESH, space = poly, matrix = 'K', order = order_dphidphi)
dphix_H1_o0, dphiy_H1_o0 = pde.h1.assemble(MESH, space = poly, matrix = 'K', order = 0)
dphix_H1_o1, dphiy_H1_o1 = pde.h1.assemble(MESH, space = poly, matrix = 'K', order = 1)
dphix_H1_order_phiphi, dphiy_H1_order_phiphi = pde.h1.assemble(MESH, space = poly, matrix = 'K', order = order_phiphi)
phi_H1b = pde.h1.assembleB(MESH, space = poly, matrix = 'M', shape = phi_H1.shape, order = order_phiphi)
phi_L2 = pde.l2.assemble(MESH, space = dxpoly, matrix = 'M', order = order_dphidphi)
D_order_dphidphi = pde.int.assemble(MESH, order = order_dphidphi)
D_order_phiphi = pde.int.assemble(MESH, order = order_phiphi)
D_order_phiphi_b = pde.int.assembleB(MESH, order = order_phiphi)
Kxx = dphix_H1 @ D_order_dphidphi @ dphix_H1.T
Kyy = dphiy_H1 @ D_order_dphidphi @ dphiy_H1.T
Cx = phi_L2 @ D_order_dphidphi @ dphix_H1.T
Cy = phi_L2 @ D_order_dphidphi @ dphiy_H1.T
D_stator_outer = pde.int.evaluateB(MESH, order = order_phiphi, edges = ind_stator_outer)
B_stator_outer = phi_H1b@ D_stator_outer @ phi_H1b.T
penalty = 1e10
J = 0; J0 = 0
for i in range(48):
J += pde.int.evaluate(MESH, order = order_phiphi, coeff = lambda x,y : j3[i], regions = np.r_[ind_trig_coils[i]]).diagonal()
# J0+= pde.int.evaluate(MESH, order = 0, coeff = lambda x,y : j3[i], regions = np.r_[ind_trig_coils[i]]).diagonal()
J = 1*J
M0 = 0; M1 = 0; M00 = 0; M10 = 0
for i in range(16):
M0 += pde.int.evaluate(MESH, order = order_phiphi, coeff = lambda x,y : m[0,i], regions = np.r_[ind_trig_magnets[i]]).diagonal()
M1 += pde.int.evaluate(MESH, order = order_phiphi, coeff = lambda x,y : m[1,i], regions = np.r_[ind_trig_magnets[i]]).diagonal()
# M00 += pde.int.evaluate(MESH, order = 0, coeff = lambda x,y : m[0,i], regions = np.r_[ind_trig_magnets[i]]).diagonal()
# M10 += pde.int.evaluate(MESH, order = 0, coeff = lambda x,y : m[1,i], regions = np.r_[ind_trig_magnets[i]]).diagonal()
aJ = phi_H1@ D_order_phiphi @J
aM = dphix_H1_order_phiphi@ D_order_phiphi @(-M1) +\
dphiy_H1_order_phiphi@ D_order_phiphi @(+M0)
aMnew = aM
print('Assembling + stuff ', time.monotonic()-tm)
Assembling + stuff 0.1913953089970164
TODO
maxIter = 100
epsangle = 1e-5
u = np.zeros(shape = Kxx.shape[0])
angleCondition = np.zeros(5)
eps_newton = 1e-8
factor_residual = 1/2
mu = 0.0001
def gss(u):
ux = dphix_H1.T@u
uy = dphiy_H1.T@u
fxx_grad_u_Kxx = dphix_H1 @ D_order_dphidphi @ sps.diags(fxx(ux,uy))@ dphix_H1.T
fyy_grad_u_Kyy = dphiy_H1 @ D_order_dphidphi @ sps.diags(fyy(ux,uy))@ dphiy_H1.T
fxy_grad_u_Kxy = dphiy_H1 @ D_order_dphidphi @ sps.diags(fxy(ux,uy))@ dphix_H1.T
fyx_grad_u_Kyx = dphix_H1 @ D_order_dphidphi @ sps.diags(fyx(ux,uy))@ dphiy_H1.T
return (fxx_grad_u_Kxx + fyy_grad_u_Kyy + fxy_grad_u_Kxy + fyx_grad_u_Kyx) + penalty*B_stator_outer
def gs(u):
ux = dphix_H1.T@u
uy = dphiy_H1.T@u
return dphix_H1 @ D_order_dphidphi @ fx(ux,uy) + dphiy_H1 @ D_order_dphidphi @ fy(ux,uy) + penalty*B_stator_outer@u - aJ + aM
def J(u):
ux = dphix_H1.T@u
uy = dphiy_H1.T@u
return np.ones(D_order_dphidphi.size)@ D_order_dphidphi @f(ux,uy) -(aJ-aM)@u + 1/2*penalty*u@B_stator_outer@u
tm = time.monotonic()
for i in range(maxIter):
gsu = gs(u)
gssu = gss(u)
w = chol(gssu).solve_A(-gsu)
# If you want to solve iteratively, use the following instead:
# ifssu = sps.diags(1/(fss(u).diagonal()), format = 'csc') #Jacobi preconditioner
# precon = lambda x : ifssu@precon_V(x)
# precon = pyamg.smoothed_aggregation_solver(fssu).aspreconditioner(cycle = 'V') #Algebraic multigrid V-cycle
# w = pde.pcg(fssu, fsu, tol = 1e-1, maxit = 100, pfuns = precon, output = True)
norm_w = np.linalg.norm(w)
norm_gsu = np.linalg.norm(gsu)
if (-(w@gsu)/(norm_w*norm_gsu) < epsangle):
angleCondition[i%5] = 1
if np.product(angleCondition)>0:
w = -gsu
print("STEP IN NEGATIVE GRADIENT DIRECTION")
else: angleCondition[i%5]=0
alpha = 1
# ResidualLineSearch
# for k in range(1000):
# if np.linalg.norm(fs(u+alpha*w)) <= np.linalg.norm(fs(u)): break
# else: alpha = alpha*factor_residual
# AmijoBacktracking
float_eps = np.finfo(float).eps
for k in range(1000):
if J(u+alpha*w)-J(u) <= alpha*mu*(-gs(u)@w) + np.abs(J(u))*float_eps: break
else: alpha = alpha*factor_residual
u = u + alpha*w
print ("NEWTON: Iteration: %2d " %(i+1)+"||obj: %2e" %J(u)+"|| ||grad||: %2e" %np.linalg.norm(fs(u))+"||alpha: %2e" % (alpha))
if(np.linalg.norm(gs(u)) < eps_newton): break
elapsed = time.monotonic()-tm
print('Solving took ', elapsed, 'seconds')
/var/folders/y9/52c8hkj57sndl5zj9s3n62n80000gn/T/ipykernel_17656/1772158371.py:4: RuntimeWarning: overflow encountered in exp
NEWTON: Iteration: 1 ||obj: -9.321500e+02|| ||grad||: 4.403518e+05||alpha: 1.250000e-01 NEWTON: Iteration: 2 ||obj: -1.230529e+03|| ||grad||: 4.136966e+05||alpha: 6.250000e-02 NEWTON: Iteration: 3 ||obj: -1.596589e+03|| ||grad||: 3.650968e+05||alpha: 1.250000e-01 NEWTON: Iteration: 4 ||obj: -1.940527e+03|| ||grad||: 2.846259e+05||alpha: 2.500000e-01 NEWTON: Iteration: 5 ||obj: -2.192195e+03|| ||grad||: 1.734370e+05||alpha: 5.000000e-01 NEWTON: Iteration: 6 ||obj: -2.243241e+03|| ||grad||: 2.063412e+05||alpha: 1.000000e+00 NEWTON: Iteration: 7 ||obj: -2.276135e+03|| ||grad||: 7.707736e+04||alpha: 1.000000e+00 NEWTON: Iteration: 8 ||obj: -2.287455e+03|| ||grad||: 2.882314e+04||alpha: 1.000000e+00 NEWTON: Iteration: 9 ||obj: -2.291545e+03|| ||grad||: 1.076473e+04||alpha: 1.000000e+00 NEWTON: Iteration: 10 ||obj: -2.292956e+03|| ||grad||: 4.010763e+03||alpha: 1.000000e+00 NEWTON: Iteration: 11 ||obj: -2.293422e+03|| ||grad||: 1.498494e+03||alpha: 1.000000e+00 NEWTON: Iteration: 12 ||obj: -2.293576e+03|| ||grad||: 5.614604e+02||alpha: 1.000000e+00 NEWTON: Iteration: 13 ||obj: -2.293621e+03|| ||grad||: 1.908430e+02||alpha: 1.000000e+00 NEWTON: Iteration: 14 ||obj: -2.293630e+03|| ||grad||: 5.017592e+01||alpha: 1.000000e+00 NEWTON: Iteration: 15 ||obj: -2.293631e+03|| ||grad||: 6.981619e+00||alpha: 1.000000e+00 NEWTON: Iteration: 16 ||obj: -2.293631e+03|| ||grad||: 1.888558e-01||alpha: 1.000000e+00 NEWTON: Iteration: 17 ||obj: -2.293631e+03|| ||grad||: 1.476924e-04||alpha: 1.000000e+00 NEWTON: Iteration: 18 ||obj: -2.293631e+03|| ||grad||: 4.214113e-10||alpha: 1.000000e+00 Solving took 1.6927945159986848 seconds
fig = MESH.pdesurf_hybrid(dict(trig = 'P1', quad = 'Q1', controls = 0), u[:MESH.np], u_height = 1)
fig.layout.scene.camera.projection.type = "orthographic"
fig.data[0].colorscale='Jet'
fig.show()